> chrismb_phy <- readRDS("/shared/statgen/microbiome/CHRISMB/processed_data/microbiome_tables/phyloseq_built.Rds")
> 
> library(chrisData)
> 
> possible_data_files <- chrisDataFile()
> 
> selected_data_files <- c(
+     "general_information",
+     "household",
+     "clinical_traits",
+     "interview",
+     "labdata",
+     "substudies",
+     "drugs_summary"
+   )
> 
> data_selected.list <- sapply(selected_data_files, function(x) chrisData(paste0(x, ".txt.gz")), USE.NAMES = T)
> 
> data_selected.joined.df <- purrr::reduce(data_selected.list, full_join, by = "AID")
> data_selected.joined.df <- data_selected.joined.df[match(chrismb_phy@sam_data$AID, data_selected.joined.df$AID),]
> 
> # Columns with all NAs must be necessarily removed 
> all_NAs <- apply(data_selected.joined.df, 2, function(x) all(is.na(x)))
> data_selected.joined.df <- data_selected.joined.df[, !all_NAs]
> 
> # Columns of class "chron" are problematic, they must be removed to be able to work with phyloseq and metadata data.frame back and forth
> data_selected.joined_NODates.df <- select(data_selected.joined.df, -contains(c("x0bp08", "x0bc03", "x0bc06")))
> 
> rm(data_selected.list)
> 
> chrismb_phy <- phy_add_metadata_variables(chrismb_phy, df = data_selected.joined_NODates.df, by = "AID")
> 
> 
> levels(chrismb_phy@sam_data$x0oh01) <- c("0", "1-9", "10-19", "20+")
> levels(chrismb_phy@sam_data$x0sm51) <- c("Never", "Former", "Current (R)", "Current (NR)")
> chrismb_phy@sam_data$x0an03a %<>% as.numeric()
> 
> chrismb_phy@sam_data$x0_ager_cat <- cut(chrismb_phy@sam_data$x0_ager, breaks = c(18, 31, 41, 51, 61, 71, max(chrismb_phy@sam_data$x0_ager)+1), include.lowest = TRUE, labels = c("18-30", "31-40", "41-50", "51-60", "61-70", "71+"))
> 
> # create some househod (hid) info
> chrismb_phy <- relocate_sample_data(chrismb_phy, x0_ager_cat, .after = x0_ager) %>% 
+   mutate_sample_data(
+     hid = paste0("h", hid),
+     hid_size_chrismb = add_N_to_catVar(hid) %>% as.character() %>% gsub("h.*\\ ", "", .) %>% parse_number(),
+   smoking_exposure_ga = 
+     case_when(
+       grepl("Current", as.character(x0sm51)) ~ "Current",
+       TRUE ~ as.character(x0sm51)
+       ) %>% 
+     factor(levels = c("Never", "Former", "Current"))
+   ) %>% 
+   relocate_sample_data(smoking_exposure_ga, .before = x0sm51)

> excluded_samples_upon_request <- readLines("/home/gantonello/CHRISMB/participants lists each study/withdrawnConsents_AIDs.txt")
> 
> chrismb_phy_core <- core(chrismb_phy, detection = threshold_detection, prevalence = threshold_prevalence)
> 
> phy_Q1 <- chrismb_phy_core %>% 
+   filter_sample_data(!(AID %in% excluded_samples_upon_request),
+                      !is.na(x0oh01),
+                      !is.na(x0_ager),
+                      !is.na(x0_sex),
+                      !is.na(x0sm51),
+                      #x0bc10 != 1, # people who smoked were reported as 1, NAs are treated as non-smoking
+                      x0dd51 == "No")
> 
> phy_Q1

> phy_Q1_genus <- phy_tax_glom(phy_Q1, "Genus")
> 
> 
> phy_Q1_meta <- meta(phy_Q1)

> dir.create("../results")

> div <- phyloseq::distance(microbiome::transform(phy_Q1_genus, "compositional"), "bray")
> 
> ord <- ordinate(physeq = phy_Q1_genus,
+                 distance = div,
+                 method = "PCoA"
+ )
> 
> fig1A <- plot_ordination(physeq = phy_Q1_genus,
+                          ordination = ord, 
+                          color = "smoking_exposure_ga")+ 
+   labs(title = NULL)+
+   theme_light()+
+   theme(legend.position = "bottom")+ 
+   scale_color_discrete(name = "Smoking habit",type = ggsci::pal_jco("default")(4)[c(1,2,4)])+ 
+   stat_ellipse(aes(group = smoking_exposure_ga), show.legend = FALSE)+
+   guides(color = guide_legend(override.aes = list(size=3)))
> 
> fig1A$labels$caption <- gsub(fig1A$labels$caption, pattern = "\n", replacement = "; ")
> fig1A$labels$x <- paste0("Principal Coordinate 1   ", strsplit(fig1A$labels$x, "   ")[[1]][[2]])
> fig1A$labels$y <- paste0("Principal Coordinate 2   ", strsplit(fig1A$labels$y, "   ")[[1]][[2]])

> dir.create("../results/Figure1/OxygenMetabolism/", recursive = T)
> 
> otutab_with_genus_name <- phy_Q1_genus %>%
+   microbiome::transform("compositional") %>%
+   abundances() %>% 
+   as.data.frame() %>% 
+   rownames_to_column("Genus") %>% # next step helps assigning the highest number of genera/families possible (because among the taxa in the table there are some families too) 
+   mutate(Genus_easier = sapply(strsplit(as.character(Genus), "_"), function(g) g[[1]]))
> 
> # get the metadata of the oxygen metabolism vs Genus
> oxygen_requirement_taxa <- read.csv("https://raw.githubusercontent.com/mcalgaro93/sc2meta/master/data/genera_methabolism.tsv", sep = "\t") %>% 
+   set_names(c("Genus_easier", "oxygen_requirements")) %>% 
+   mutate(oxygen_requirements = gsub("F ", "Facultative ", oxygen_requirements, fixed = T))
> 
> # merge otu table and oxygen requirement data by simplified genus name
> molten_df_merged_oxygen <- merge(otutab_with_genus_name, oxygen_requirement_taxa, by = "Genus_easier", all.x = TRUE)
> molten_df_merged_oxygen$oxygen_requirements <- ifelse(is.na(molten_df_merged_oxygen$oxygen_requirements), "NA",molten_df_merged_oxygen$oxygen_requirements)
> 
> # aggregate by oxygen metabolism, instead of genera, to make code easier
> 
> tmp <- molten_df_merged_oxygen %>% 
+   dplyr::select(!contains("Genus")) %>%
+   group_by(oxygen_requirements) %>% 
+   summarise_all(sum) %>%
+   column_to_rownames("oxygen_requirements") %>% 
+   t() %>% as.data.frame() %>% 
+   rownames_to_column("sample_name")
> 
> oxygen_molten_metadata <- tmp %>% reshape2::melt() %>% 
+   merge(., phy_Q1_meta, all.x = TRUE, by = "sample_name") %>% 
+   dplyr::rename(oxygen = variable)
> 
> saveRDS(oxygen_molten_metadata, "../results/Figure1/OxygenMetabolism/oxygen_with_metadata_df.RDS")

> library(rstatix)
> 
> signif <- filter(oxygen_molten_metadata, oxygen != "NA") %>%
+   mutate(aerobiosis = oxygen) %>% 
+   group_by(oxygen) %>% 
+   rstatix::wilcox_test(value ~ smoking_exposure_ga) %>% 
+   adjust_pvalue(method = "BH") %>% 
+   #add_significance("p.adj") %>% 
+   add_xy_position(x = "oxygen",step.increase = 0.05) %>% 
+   filter(p.adj < 0.05)
> 
> fig1C <- ggboxplot(data = filter(oxygen_molten_metadata, oxygen != "NA"),
+                    x = "oxygen", 
+                    y = "value", 
+                    fill = "smoking_exposure_ga") + 
+   ylab("Relative Abundance") + 
+   xlab("")+ 
+   scale_fill_discrete(name = "Smoking habit", type = ggsci::pal_jco()(4)[c(1,2,4)]) + 
+   stat_pvalue_manual(signif, tip.length = 0.01) + 
+   theme_light()+
+   theme(legend.position = "none")

> if(!file.exists("../results/Figure2/deseq_smoking_regression_datafames.rds")){
+ deseq_raw_results_Q2 <- list()
+ for(gr in c("x0sm61_bin5cont", "x0sm61_bin2cont", "x0sm61")){
+   
+   # create deseq object
+   fla <- as.formula(paste0("~ x0_ager + x0_sex + x0oh01 + ", gr))
+   
+   phy_q2_genus_deseq_q2 <- phyloseq_to_deseq2(phy_Q2_genus, design = fla)
+   
+   # calculate geometric means, important to center the counts
+   
+   geom_means_q2 <- apply(counts(phy_q2_genus_deseq_q2), 1, function(x)
+       exp(mean(log(x[x > 0]))))
+   
+   # estimate size factors, important to account for library size variations
+   estim_size_q2 <- estimateSizeFactors(phy_q2_genus_deseq_q2, geoMeans = geom_means_q2)
+   
+   # compute the DESeq2 glm, accounting for the estimate size factor and geometric means
+   
+   deseq2_results_q2 <- DESeq(estim_size_q2, fitType="local")
+   deseq_raw_results_Q2[[gr]] <- deseq2_results_q2
+ }
+ 
+ taxonomy <- as.data.frame(tax_table(phy_Q2_genus))[c("Genus", "Phylum")]
+ q2_results <- lapply(deseq_raw_results_Q2, function(res) results(res) %>% 
+                        as.data.frame() %>% 
+                        rownames_to_column("Genus") %>%
+                        merge(taxonomy, by = "Genus", sort = FALSE, all.x = TRUE)
+ )
+ 
+ saveRDS(q2_results, "../results/Figure2/deseq_smoking_regression_datafames.rds")
+ rm(taxonomy)
+ 
+ } else{
+   
+   q2_results <- readRDS("../results/Figure2/deseq_smoking_regression_datafames.rds")
+ }

> Fig2BCD <-
+   filter(
+     oxygen_molten_metadata,
+     oxygen != "NA" &
+       smoking_exposure_ga == "Current",
+     x0sm61 < 60
+   ) %>%
+   group_split(oxygen) %>%
+   lapply(
+     function(ox)
+       ggscatter(ox,
+                 x = "x0sm61",
+                 y = "value",
+                 alpha = 0.5,
+                 color = "gray30") +
+       geom_smooth(se = TRUE) +
+       theme_light()+
+       labs(
+         x = "Tobacco smoked per day (g)",
+         y = "Relative Abundance",
+         title = unique(ox$oxygen) 
+         )
+ )

> filter(
+     oxygen_molten_metadata,
+     oxygen == "Aerobic" &
+       smoking_exposure_ga == "Current",
+     x0sm61 < 60
+   ) %>% 
+   lm(value ~ I(1/x0sm61) + x0_ager + x0_sex + x0oh01, data = .) %>% 
+   summary() %>% 
+   broom::tidy() %>%
+   mutate(term = c("Intercept", "1/tobacco exposure intensity (g/day)", "Age", "Sex", "1-9 teeth", "10-19 teeth", "20+ teeth")) %>% 
+   kbl(caption = "The relative abundance of AEROBIC bacteria decreases with increasing grams of tobacco smoked in a day, up to 10 g/day. After 10 grams, it is no longer statistically associated (data not shown)", digits = 3) %>% 
+   kable_styling()

> library(DESeq2)
> 
> phy_q3_genus_deseq <- phyloseq_to_deseq2(phy_Q3_genus, design = ~ x0_ager_cat + x0_sex + x0oh01 + years_since_quitting_or_reducing)
> 
> # calculate geometric means, important to center the counts
> geom_means_q3 <- apply(counts(phy_q3_genus_deseq), 1, function(x)
+       exp(mean(log(x[x > 0]))))
> 
> # estimate size factors, important to account for library size variations
> estim_size_q3 <- estimateSizeFactors(phy_q3_genus_deseq, geoMeans = geom_means_q3)
> 
> # compute the DESeq2 glm, accounting for the estimate size factor and geometric means
> 
> deseq2_results_q3 <- DESeq(estim_size_q3, fitType="local")
> 
> results_q3 <- results(deseq2_results_q3, name = "years_since_quitting_or_reducing")
> 
> results_q3.df <- as.data.frame(results_q3)
> 
> DT::datatable(arrange(results_q3.df, padj) %>% mutate_if(is.numeric, round, 3), caption = "Table 'Years since quitting' | Results of regression of microbial genera in response to the years since quitting smoking. The model was corrected for age, sex and number of teeth, nominal P-values correceted for False Discovery Rate with the Benjamini-Hochberg method.")

> fig3BCD <-
+   filter(
+     oxygen_molten_metadata,
+     oxygen != "NA" &
+       sample_name %in% sample_names(phy_Q3_genus) &
+     years_since_quitting_or_reducing <= 40
+   ) %>%
+   group_split(oxygen) %>%
+   lapply(
+     function(ox)
+       ggscatter(ox,
+                 x = "years_since_quitting_or_reducing",
+                 y = "value",
+                 alpha = 0.5
+       ) +
+       geom_smooth(se = TRUE) +
+       xlab("Years since quitting smoking") +
+       ylab("Relative Abundance") +
+       ggtitle(unique(ox$oxygen)) +
+       
+       theme_light()+geom_point(alpha = 0.1, color = "gray50")
+   )

> filter(
+     oxygen_molten_metadata,
+     oxygen == "Aerobic" &
+         sample_name %in% sample_names(phy_Q3_genus) &
+         years_since_quitting_or_reducing <= 20) %>% 
+     lm(value ~ years_since_quitting_or_reducing + x0_ager + x0_sex + x0oh01, data = .) %>% 
+   summary() %>% 
+   broom::tidy() %>%
+   mutate(term = c("Intercept", "Years since quitting", "Age", "Sex", "1-9 teeth", "10-19 teeth", "20+ teeth")) %>% 
+   kbl(caption = "The relative abundance of AEROBIC bacteria increases with increasing years since quitting, up to 20 years. After 20 years, it is no longer statistically associated (data not shown)", digits = 3) %>% 
+   kable_styling()

> filter(
+     oxygen_molten_metadata,
+     oxygen == "Anaerobic" &
+         sample_name %in% sample_names(phy_Q3_genus) &
+         years_since_quitting_or_reducing <= 20) %>% 
+     lm(value ~ years_since_quitting_or_reducing + x0_ager + x0_sex + x0oh01, data = .) %>% 
+   summary() %>% 
+   broom::tidy() %>%
+   mutate(term = c("Intercept", "Years since quitting", "Age", "Sex", "1-9 teeth", "10-19 teeth", "20+ teeth")) %>% 
+   kbl(caption = "The relative abundance of ANAEROBIC bacteria increases with increasing years since quitting, up to 20 years. After 20 years, it is no longer statistically associated (data not shown)", digits = 3) %>% 
+   kable_styling()

> dir.create("../results/Figure4/DiffAbundGenera/DESeq2/", recursive = T)
> 
> if(!file.exists("../results/Figure4/DiffAbundGenera/DESeq2/figure4_deseq2_dataframe.list.RDS")){
+   
+ library(DESeq2)
+ #############################################################################
+ 
+ deseq_object_picrust <- phyloseq(otu_table(picrust_pathways_matrix_Q1, taxa_are_rows = TRUE), sample_data(phy_Q1))  %>% 
+   core(detection = 10, prevalence = 0.05) %>% 
+   phyloseq_to_deseq2(design = ~ x0_sex + x0_ager_cat + x0oh01 + smoking_exposure_ga)
+ 
+ geom_means <- apply(counts(deseq_object_picrust), 1, function(x) exp(mean(log(x[x>0]))))
+ 
+ # estimate size factors, important to account for library size variations
+ estim_size <- estimateSizeFactors(deseq_object_picrust, geoMeans = geom_means)
+ 
+ # compute the DESeq2 glm, accounting for the estimate size factor and geometric means
+ 
+ deseq2_results_picrust <- DESeq(estim_size, fitType = "local")
+ 
+ saveRDS(deseq2_results_picrust, "../results/Figure4/DiffAbundGenera/DESeq2/figure4_deseq2_dataframe.list.RDS")
+ } else{
+   deseq2_results_picrust <- readRDS("../results/Figure4/DiffAbundGenera/DESeq2/figure4_deseq2_dataframe.list.RDS")
+ }

> dir.create("../results/Figure4/DiffAbundGenera/ALDEx2/", recursive = T)
> 
> if (!file.exists("../results/Figure4/DiffAbundGenera/ALDEx2/ALDEx2.DA.rawResults.list.RDS")) {
+   
+   library(ALDEx2)
+   
+   model_mtx <-
+     model.matrix(
+       ~  x0_sex + x0_ager_cat + x0oh01 + smoking_exposure_ga,
+       data = dplyr::select(
+         phy_Q1_meta,
+         x0oh01,
+         x0_sex,
+         x0_ager_cat,
+         smoking_exposure_ga
+       ) %>% .[intersect(colnames(picrust_pathways_matrix_Q1),
+                         rownames(phy_Q1_meta)), ]
+     )
+   
+   aldex.glm.pathway_smoking_ga <-
+     picrust_pathways_matrix_Q1[, intersect(colnames(picrust_pathways_matrix_Q1),
+                                                      rownames(phy_Q1_meta))] %>% 
+     aldex.clr(
+               model_mtx,
+               mc.samples = 150,
+               denom = "all", 
+               useMC = T
+               )
+   
+   glm.test_pathway_smoking_ga <-
+     aldex.glm(aldex.glm.pathway_smoking_ga, model_mtx)
+   
+   glm.effects_pathway_smoking_ga <-
+     aldex.glm.effect(aldex.glm.pathway_smoking_ga, verbose = FALSE, include.sample.summary = T, useMC = T)
+   
+   aldex.results <- list(clr.model = aldex.glm.pathway_smoking_ga,
+                         glm.test = glm.test_pathway_smoking_ga,
+                         eff_size.test = glm.effects_pathway_smoking_ga)
+   
+   
+   saveRDS(aldex.results, "../results/Figure4/DiffAbundGenera/ALDEx2/ALDEx2.DA.rawResults.list.RDS")
+ 
+   detach("package:ALDEx2", unload = T)
+   detach("package:zCompositions", unload = T)
+   
+ } else{
+   aldex.results <- readRDS("../results/Figure4/DiffAbundGenera/ALDEx2/ALDEx2.DA.rawResults.list.RDS")
+ }

> dir.create("../results/Figure4/DiffAbundGenera/ANCOM-BC/", recursive = T)
> 
> if (!file.exists("../results/Figure4/DiffAbundGenera/ANCOM-BC/ANCOM-BC.DA.rawResults.list.RDS")) {
+ tmp <- phyloseq(otu_table(picrust_pathways_matrix_Q1, taxa_are_rows = T), sample_data(phy_Q1_genus %>%
+   select_sample_data(x0_sex, x0_ager_cat, x0oh01, smoking_exposure_ga) %>% meta())) %>%
+   select_sample_data(x0_sex, x0_ager_cat, x0oh01, smoking_exposure_ga) %>% 
+   core(detection = 10, prevalence = 0.05) %>%
+   mia::makeTreeSummarizedExperimentFromPhyloseq()
+ 
+ 
+ library(ANCOMBC)
+ 
+ ancom_results_picrust <- ancombc2(data = tmp,
+         assay_name = "counts", 
+         tax_level = NULL,
+         fix_formula = "x0_sex + x0_ager_cat + x0oh01 + smoking_exposure_ga", 
+         group = "smoking_exposure_ga",
+         struc_zero = TRUE,
+         iter_control = list(tol = 1e-5, max_iter = 20, verbose = FALSE),
+         em_control = list(tol = 1e-5, max_iter = 20),  
+         p_adj_method = "BH",
+         n_cl = 1,
+         alpha = 0.05)
+ 
+ saveRDS(ancom_results_picrust, "../results/Figure4/DiffAbundGenera/ANCOM-BC/ANCOM-BC.DA.rawResults.list.RDS")
+ 
+ }else{
+   
+   ancom_results_picrust <- readRDS("../results/Figure4/DiffAbundGenera/ANCOM-BC/ANCOM-BC.DA.rawResults.list.RDS")
+   
+ }